﻿
# y = a^(-1) . z
# y -> [4x1],	a -> [4x4],	z -> [4x1]



a = [
		[rand, rand, rand, rand],
		[rand, rand, rand, rand],
		[rand, rand, rand, rand],
		[rand, rand, rand, rand]
	  ]

z = [ rand, rand, rand, rand ]
	  
def lu_invert(a, z)
	t1 = a[0][0]*a[1][1]-a[0][1]*a[1][0]
	t2 = a[0][0]*a[1][3]-a[0][3]*a[1][0]
	t3 = a[2][2]*a[3][3]-a[2][3]*a[3][2] 
	t4 = a[2][1]*a[3][2]-a[2][2]*a[3][1]
	t5 = a[0][2]*a[1][0]-a[0][0]*a[1][2] 
	t6 = a[0][2]*a[1][3]-a[0][3]*a[1][2] 
	t7 = a[0][2]*a[1][1]-a[0][1]*a[1][2] 
	t8 = a[0][1]*a[1][3]-a[0][3]*a[1][1] 
	t9 = a[2][0]*a[3][3]-a[2][3]*a[3][0] 
	t10 = a[2][0]*a[3][2]-a[2][2]*a[3][0] 
	t11 = a[2][1]*a[3][0]-a[2][0]*a[3][1] 
	t12 = a[2][3]*a[3][1]-a[2][1]*a[3][3] 
	t13 = a[2][0]*z[3]-a[3][0]*z[2]
	t14 = a[0][3]*a[2][0]-a[0][0]*a[2][3] 
	t15 = a[1][0]*a[2][3]-a[1][3]*a[2][0] 
	t16 = a[2][2]*z[3]-a[3][2]*z[2] 
	t17 = a[2][3]*z[3]-a[3][3]*z[2] 
	t18 = a[3][1]*z[2]-a[2][1]*z[3] 
	t19 = a[0][2]*a[2][1]-a[0][1]*a[2][2] 
	t20 = a[0][1]*a[2][0]-a[0][0]*a[2][1] 
	t21 = a[0][1]*a[2][3]-a[0][3]*a[2][1] 
	t22 = a[0][3]*a[2][2]-a[0][2]*a[2][3] 
	t23 = a[0][0]*a[2][2]-a[0][2]*a[2][0] 
	t24 = a[1][2]*a[2][3]-a[1][3]*a[2][2] 
	t25 = a[1][1]*a[2][2]-a[1][2]*a[2][1] 
	t26 = a[1][3]*a[2][1]-a[1][1]*a[2][3] 
	t27 = a[1][0]*a[2][1]-a[1][1]*a[2][0] 
	t28 = a[1][2]*a[2][0]-a[1][0]*a[2][2] 
	t29 = a[3][2]*z[0] 
	t30 = a[3][3]*z[0] 
	t31 = a[3][1]*z[0] 
	t32 = a[3][0]*z[0] 
	t33 = a[3][1]*z[1] 
	t34 = a[3][2]*z[1] 
	t35 = a[3][3]*z[1] 
	t36 = a[3][0]*z[1] 

	d = t7*t9+t10*t8+t11*t6+t12*t5-t2*t4-t1*t3 

	[
		-(t16*t8+t17*t7+t18*t6+t19*t35+t21*t34+t22*t33+t24*t31+t25*t30+t26*t29)/d,
		-(t13*t6-t17*t5-t22*t36+t23*t35+t14*t34-t24*t32+t28*t30+t15*t29-t16*t2)/d,
		(t13*t8+t21*t36-t20*t35+t14*t33+t26*t32+t15*t31-t27*t30+t18*t2+t1*t17)/d,
		(t13*t7+t18*t5+t19*t36+t20*t34+t23*t33+t25*t32+t28*t31+t27*t29-t1*t16)/d
	]
end

y = lu_invert a,z
puts "A: #{a}\n\n", "Z: #{z}\n\n", "A.inv * Z = #{y}\n\n"
